#define _USE_MATH_DEFINES
#include <cmath>
#include "fluid.h"
#include "pipe.h"
#include "solver.h"
#include <iostream>
#include <fstream>
#include <iomanip>

const double PI = 3.14159265358979323846;

void printResults(const MultiphaseFlowSolver::Results& results) {
    std::cout << std::fixed << std::setprecision(3);
    std::cout << "\nResults along the pipe:\n";
    std::cout << "Position(m) Pressure(MPa) Temp(C) GasHold OilHold WaterHold Vg(m/s) Vo(m/s) Vw(m/s)\n";
    
    int numPoints = results.pressure.size();
    double dx = 4000.0 / (numPoints - 1);
    
    // 打印几个代表性位置的结果
    for (int i = 0; i < numPoints; i += numPoints/10) {  // 每10%打印一次
        double x = i * dx;
        double tempC = results.temperature[i] - 273.15;  // 转换为摄氏度
        
        std::cout << std::setw(10) << x << " "
                  << std::setw(10) << results.pressure[i]/1e6 << " "  // 转换为MPa
                  << std::setw(8) << tempC << " "  // 已经转换为摄氏度
                  << std::setw(7) << results.gasHoldup[i] << " "
                  << std::setw(7) << results.oilHoldup[i] << " "
                  << std::setw(9) << results.waterHoldup[i] << " "
                  << std::setw(7) << results.gasVelocity[i] << " "
                  << std::setw(7) << results.oilVelocity[i] << " "
                  << std::setw(7) << results.waterVelocity[i] << "\n";
    }
}
int main() {
    try {
        std::cout << "Creating pipe parameters..." << std::endl;
        PipeParameters params;
        // 根据表格数据设置参数
        params.length = 4000.0;          // 管道长度(m)
        params.diameter = 0.12;          // 管道内径(m) 
        params.wallThickness = 0.049;    // 管壁厚度(m)
        params.roughness = 8.0e-5;       // 管道粗糙度(m)
        params.massFlowRate = 12.0;      // 管道输量(kg/s)
        params.heatTransferCoef = 1.0;   // 总传热系数(W/(m2·K))
        
        // 压力和温度参数
        params.inletPressure = 5.0e6;    // 入口压力5MPa
        params.outletPressure = 5.0e5;   // 出口压力0.5MPa
        params.inletTemp = 323.15;       // 入口温度50℃
        params.ambientTemp = 279.15;     // 环境温度6℃
        
        // 高程参数
        params.startElevation = 0.0;     // 起点高程(m)
        params.endElevation = 0.0;       // 终点高程(m)
        
        // 相分数和流速
        params.inletGasVolumeFraction = 0.3;  // 入口气相体积分数
        
        // 从PVT数据中获取初始密度
        double rho_g = 1.15481;          // 标准状况下气相密度(kg/m3)
        double rho_l = 793.895;          // 标准状况下油相密度(kg/m3)
        double rho_w = 999.05;           // 标准状况下水相密度(kg/m3)
        
        // 计算截面积
        double area = PI * params.diameter * params.diameter / 4.0;
        
        // 根据总质量流量计算各相流量
        double total_mass_flow = params.massFlowRate;
        double water_cut = 0.05;         // 含水率5%
        
        // 计算各相质量流量
        double gas_mass_flow = total_mass_flow * params.inletGasVolumeFraction * rho_g / 
            (params.inletGasVolumeFraction * rho_g + (1-params.inletGasVolumeFraction) * 
            ((1-water_cut)*rho_l + water_cut*rho_w));
            
        double liquid_mass_flow = total_mass_flow - gas_mass_flow;
        double water_mass_flow = liquid_mass_flow * water_cut;
        double oil_mass_flow = liquid_mass_flow * (1 - water_cut);
        
        // 计算各相速度
        params.inletGasVelocity = gas_mass_flow / (area * params.inletGasVolumeFraction * rho_g);
        params.inletLiquidVelocity = liquid_mass_flow / 
            (area * (1-params.inletGasVolumeFraction) * ((1-water_cut)*rho_l + water_cut*rho_w));
        
        // 设置时间步长参数（根据表格数据）
        params.maxTimeStep = 5.0;        // 最大模拟时间步长(s)
        params.minTimeStep = 0.01;       // 最小模拟时间步长(s)
        params.simulationTime = 60.0;  // 模拟时间(s)

        std::cout << "Creating fluid..." << std::endl;
        std::string pvtFile = "../test_case/test3.tab";
        
        // 检查文件是否存在
        std::ifstream f(pvtFile.c_str());
        if (!f.good()) {
            std::cerr << "Error: Cannot find file: " << pvtFile << std::endl;
            return 1;
        }

        Fluid fluid(pvtFile);
        
        // 在创建solver之前添加参数检查
        if (params.inletPressure <= 0 || params.outletPressure <= 0) {
            std::cerr << "Error: Invalid pressure values" << std::endl;
            return 1;
        }

        if (params.inletTemp <= 0 || params.ambientTemp <= 0) {
            std::cerr << "Error: Invalid temperature values" << std::endl;
            return 1;
        }

        if (params.inletGasVolumeFraction < 0 || params.inletGasVolumeFraction > 1) {
            std::cerr << "Error: Invalid gas volume fraction" << std::endl;
            return 1;
        }

        std::cout << "Creating solver..." << std::endl;
        std::ofstream detailedFile("detailed_results.txt");
        if (detailedFile.is_open()) {
            detailedFile << "Time(s) Length(m) Pressure(Pa) Temperature(K) GasHoldup(-) "
                         << "LiquidHoldup(-) WaterHoldup(-) GasVelocity(m/s) "
                         << "LiquidVelocity(m/s) WaterVelocity(m/s) GasDensity(kg/m3) "
                         << "LiquidDensity(kg/m3) WaterDensity(kg/m3) GasViscosity(Pa.s) "
                         << "LiquidViscosity(Pa.s) WaterViscosity(Pa.s)\n";
            detailedFile.close();
        }
        MultiphaseFlowSolver solver(fluid, params);
        
        std::cout << "Starting solution..." << std::endl;
        auto results = solver.solve();
        
        std::cout << "Solution completed successfully." << std::endl;
        printResults(results);
        
        // 保存结果到文件
        std::ofstream outFile("results.txt");
        if (outFile.is_open()) {
            outFile << "x(m),P(Pa),T(K),alpha_g,alpha_o,alpha_w,v_g(m/s),v_o(m/s),v_w(m/s)\n";
            for (size_t i = 0; i < results.pressure.size(); i++) {
                outFile << i * params.length / (results.pressure.size()-1) << ","
                       << results.pressure[i] << ","
                       << results.temperature[i] << ","
                       << results.gasHoldup[i] << ","
                       << results.oilHoldup[i] << ","
                       << results.waterHoldup[i] << ","
                       << results.gasVelocity[i] << ","
                       << results.oilVelocity[i] << ","
                       << results.waterVelocity[i] << "\n";
            }
            outFile.close();
            std::cout << "\nDetailed results have been saved to 'detailed_results.txt'\n";
        }
        return 0;
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << std::endl;
        return 1;
    }
}